function CC_IV(spike_flag)
% Analyze a CCIV protocol
% Modified (heavily) from POK version, 11/26 and 11/27/98 P. Manis
% This returns the modified control worksheet
% since some analysis results are placed in the sheet, and the data are all
% corrected in advance (this used to call prepare data...)
% added computation of inverted boltzman for analysis of HWF data. 1/7/99 P. Manis
% plot_flag: 0 means plot only, 1 means print, 2 means print and close window
% negative means no display.

% n is flag for spike analysis - omit if n is 1, do if n is 2

% add return of ivresult structure 3/2/99. ivresult contains:
% ivresult.timestamp = time of latest analysis
% ivresult.iss is mean current
% ivresult.vss is mean voltage
% ivresult.iss_n is mean current excluding spikes
% ivresult.vss_n is mean voltage excluding traces with spikes
% ivresult.vpk_n is peak voltage excluding traces with spikes
%
% Modified 2/10/2000 : changed so that changing rates in the middle of the
% run does not cause subsequent traces to be incorrectly analyzed.
% i.e., all time pointers are recomputed for every trace, depending on the
% current rate...
% Affected also: find_spikes, spike_sum, etc.
%
% added test for hyperpolarizing voltages: measure b/a ratio (ss/peak), and
% the tau decay from the peak of the hyperpolarizaiotn...
% This allows us to distinguish T from D stellate cells in the VCN.
% the result is printed in the window, and stored in the database.

%try
sf = getmainselection;
if(sf > 0)
    pflag = getplotflag;
    QueMessage('CC analysis', 1); % clear the que
    for i = 1:length(sf)
        cc_iv2(sf(i), pflag, spike_flag);
    end;
end;
%catch
%  watchoff;
%   QueMessage('Error in Current-clamp Analysis routine', 1);
%end;


function cc_iv2(sf, plot_flag, n)

global ALLCH DFILE VOLTAGE CURRENT
global CONTROL

err = 0; %#ok<NASGU>
fsize = 7; %#ok<NASGU>
msize = 3; %#ok<NASGU>

QueMessage('CCIV - Starting');
dat = [];
time = [];
ivresult = []; % initialize it
do_spike = 1;
if(nargin == 3)
    do_spike = number_arg(n)-1;	% 1 turns off detailed spike analysis; 2 turns it on
else
    do_spike = 0; % assume we do not want it
end;

[DFILE, err] = analysis_setup(DFILE, sf); % get data and parameters.

if(err ~= 0)
    return;
end;

vch = 1;
ich = 2;
if(DFILE.nr_channel(1) <= 3)
    %   VOLTAGE = ALLCH{vch}-number_arg(CONTROL(sf).jpot);  %% note junctionj
    %   potential is added in scaledat, which always happens before we get
    %   here, to VOLTAGE.
    %  CURRENT = ALLCH{ich};
else
    % ask which cell you want (and the channel assignments to use!)
    chansel = inputdlg({'V CH', 'I CH'}, 'CC_IV Channel Selection for MultiChannel Recordings', 1, ...
        {num2str(3), num2str(4)});
    if(isempty(chansel))
        return; % probably canceled.
    end;
    vch = str2num(chansel{1});
    ich = str2num(chansel{2});
    VOLTAGE = ALLCH{vch}-number_arg(CONTROL(sf).jpot); % jpot is not added in ALLCH.
    CURRENT = ALLCH{ich};
end;


% abstract general information
spike_thresh=number_arg(CONTROL(sf).thresh);
QueMessage(sprintf('CCIV: Spike Threshold: %7.2f mV', spike_thresh));
protocol=deblank(lower(CONTROL(sf).protocol));

[records,pts]=size(CURRENT);

% get windows
ts1=number_arg(CONTROL(sf).durho);
ts2=number_arg(CONTROL(sf).durs1);
ts3=number_arg(CONTROL(sf).durs2);
ts4=number_arg(CONTROL(sf).durs3);
tsrmp = ts1;
hypflag = 0;
all_ivs = {'cciv', 'cciv-l', 'ccivl', 'ap-iv', 'ap_iv', 'ap-iv2', 'ap-ivl', 'ap-ivl2', 'drugtestiv'};
all_hyps = {'hyp', 'hyps', 'hypsave', 'hypl', 'hypwatch', 'ap-hyp', 'ap_hyp', 'ap_hypsynch', 'ap-hyp2', ...
    'ap-hyp3', 'mc-hyp2', 'dc-hyp2'};
% designate times according to protocol
switch(protocol)
    case all_ivs % standard IV from -60 mV
        tsrmp = ts1;
        tdel=ts1; tstepdur=ts2;
        tsl = ts1; % spike latency reference point
        tsw1 = ts1;
        tsw2 = ts2;
    case 'iv'
        tsrmp = ts1;
        tdel=ts1; tstepdur=ts2;
        tsl = ts1;
        tsw1 = ts1;
        tsw2 = ts2;
    case 'ivf'
        tsrmp = ts1;
        tdel=ts1; tstepdur=ts2;
        tsl = ts1;
        tsw1 = ts1;
        tsw2 = ts2;
    case 'drugtestiv'
        tsrmp = 10;
        tdel = 210;
        tstepdur = 500;
        ts1 = tdel;
        tsw1 = tdel;
        tsw2 = tstepdur+tdel;
    case all_hyps % variable prepulse level before evoking spikes
        tsrmp = ts1;
        tstepdur=ts2; % measure IV during prepulse
        tsw1 = ts1+ts2;
        tdel = ts1;
        tsw2 = ts1+ts2+ts3; % but do spike analysis/counts in second window
        tsl = tsw1;
        hypflag = 1; % set so we can do boltzman fit...
    otherwise
        QueMessage(sprintf('CC_IV - Error: protocol %s not recoginized - Exiting', protocol));
        return;
end

% compute the time base for plotting (time is [rec, npts] array for EACH record)
time=make_time(DFILE);
tmax=max(max(time));
RATES = (DFILE.rate .* DFILE.nr_channel) / 1000; % array of sampling rates, convert to msec

% mean times:
m1=tdel+tstepdur*0.5;
m2=tdel+tstepdur*0.99;
m3=tdel+tstepdur*0.97; % use for h current assessment

%peak times:
p1=tdel+tstepdur*0.01;
p2=tdel+tstepdur*0.50;
p3=tdel+tstepdur*0.25;


% calculate measurement times
trmp = floor(tsrmp./RATES); % for RMP/ihold determination

mean_t1=floor(m1./RATES); % for ss iv
mean_t2=floor(m2./RATES);
mean_t3=floor(m3./RATES);

peak_t1=floor(p1./RATES); % for peak iv
peak_t2=floor(p2./RATES);
peak_tfit=floor(p3./RATES);

sw1 = floor(tsw1./RATES);
sw2 = floor(tsw2./RATES);


QueMessage('CCIV - Finding mean and peak voltage');
% measure holding current and "rmp"
hold_cur=zeros(records,1);
rmp_volt=hold_cur;
mean_volt=hold_cur;
final_volt=hold_cur;
mean_cur = hold_cur;
for i=1:records
    hold_cur(i) = mean(CURRENT(i,1:trmp(i))); % save in array for later usage
    rmp_volt(i) = mean(VOLTAGE(i,1:trmp(i)));

    % measure steady-state IV and voltage
    mean_volt(i) = mean(VOLTAGE(i,mean_t1(i):mean_t2(i)),2);
    final_volt(i) = mean(VOLTAGE(i,mean_t3(i):mean_t2(i)),2);
    mean_cur(i) = mean(CURRENT(i,mean_t1(i):mean_t2(i)),2);
end

CONTROL(sf).iHold = mean(hold_cur);
CONTROL(sf).Rmp = mean(mean(VOLTAGE(:,1:trmp)));

% now sort the arrays by mean current... important This must be carrie through later (first spike, etc.)
[mean_cur, mcl] = sort(mean_cur);
final_volt=final_volt(mcl);
mean_volt=mean_volt(mcl);

% but first smooth the voltage out a bit
vsmo=zeros(records, size(VOLTAGE, 2));
for i = 1:records
    fsamp = 1000/RATES(i); % get sampling frequency
    fco = 3000;		% cutoff frequency in Hz
    wco = fco/(fsamp/2); % wco of 1 is for half of the sample rate, so set it like this...
    if(wco < 1) % if wco is > 1 then this is not a filter!
        [b, a] = butter(8, wco); % butterworth
        vsmo(i,:) = filter(b, a, VOLTAGE(i,:)); % filter all the traces...
    else
        vsmo(i,:) = VOLTAGE(i,:);
    end
end

% find peak IV voltage (current is same as steady-state)
peak_volt=zeros(records, 1);
peak_v_index=peak_volt;
peak_tlist=peak_volt;
for i = 1:records
    [peak_volt(i), peak_v_index(i)] = min(vsmo(i,peak_t1(i):peak_t2(i)));
    peak_v_index(i) = peak_v_index(i) + peak_t1(i);
    peak_tlist(i) = peak_v_index(i)*RATES(i);
end
peak_volt=peak_volt(mcl);
peak_tlist=peak_tlist(mcl); % sort by current

% get spike times
QueMessage('CCIV - Finding spikes');
[first_spike, first_isi, nr_spikes, spike_train]=find_spikes(DFILE, tsw1, tsw2, spike_thresh);
% disp(sprintf('Rec: found %4d spikes', sum(nr_spikes)))
first_isi=max(first_isi,0);
%
% sort first spike list by mean current list...
first_spike=first_spike(mcl);
first_isi=first_isi(mcl);
nr_spikes=nr_spikes(mcl);

% make lists that include only spikes during the stimulus window
% for spike times:
mcurfsl = mean_cur(first_spike>0);
mvoltfsl = mean_volt(first_spike>0);
[mvoltfsl, mvs] = sort(mvoltfsl);
fsl = first_spike(first_spike>0);
fsl = fsl(mvs);
mcurfisi = mean_cur(first_isi>0);
mvoltfisi = mean_volt(first_isi>0);
fisi = first_isi(first_isi>0);
v = size(spike_train);
ISILAT=[]; % initialize the ISILAT structure
for i=1:v(2)
    ISILAT(i).nspikes = length(spike_train(i).latency);
    ISILAT(i).isi= diff(spike_train(i).latency);
    ISILAT(i).lat = spike_train(i).latency;
    ISILAT(i).source = spike_train(i).source; % save the source trace number too
end

deadwin = CONTROL(sf).deadwin; % dead window for most analyses EXCEPT spikes, which have their own values.

% analyze the spikes - hw, dv/dt, amplitude from rest, etc.
QueMessage('CCIV - Analyzing Spikes');
% initialize the structure
spike.timestamp = datestr(now);
spike.edit = 0; % not editable...
spike.n = 0;
spike.measthresh = NaN; % mean(spikesumm.measthresh); % fraction of height from rest to spike peak where measurement was made
spike.vthresh = NaN; % mean(spikesumm.vthr); % true threshold (2nd derivative method)
spike.vthrmin = NaN; % min(spikesumm.vthr); % min value (first spike...)
spike.deadwin = NaN; % mean(spikesumm.deadwin);
spike.minisi = NaN; % mean(spikesumm.minisi);
spike.max= NaN; % mean(spikesumm.vmax); spike.max_sd = std(spikesumm.vmax);
spike.apht = NaN; % mean(spikesumm.apht); spike.apht_sd = std(spikesumm.apht);
spike.half_ht = NaN; % mean(spikesumm.half_ht); spike.half_ht_sd = std(spikesumm.half_ht);
spike.min=NaN; % mean(spikesumm.vmin); spike.min_sd = std(spikesumm.vmin);
spike.dvmax = NaN; % mean(spikesumm.dvmax); spike.dvmax_sd = std(spikesumm.dvmax);
spike.dvmin = NaN; % mean(spikesumm.dvmin); spike.dvmin_sd = std(spikesumm.dvmin);
%      dvvmax = spikesumm.dvvmax(find(~isnan(spikesumm.dvvmax)));
spike.dvvmax = NaN; % mean(dvvmax); spike.dvvmax_sd = std(dvvmax);
%      dvvmin = spikesumm.dvvmin(find(~isnan(spikesumm.dvvmin)));
spike.dvvmin = NaN; % mean(dvvmin); spike.dvvmin_sd = std(dvvmin);
%     hwa = spikesumm.hwa(find(~isnan(spikesumm.hwa))); hwb = spikesumm.hwb(find(~isnan(spikesumm.hwb)));
spike.hwa = NaN; % mean(hwa); spike.hwa_sd = std(hwa); % eliminate the ones where coulnd't get a measurement
spike.hwb = NaN; % mean(hwb); spike.hwb_sd = std(hwb);
spike.ahpdepth =  NaN; % mean(spikesumm.ahpdepth); spike.ahpdepth_sd = std(spikesumm.ahpdepth);

CONTROL(sf).spike=spike; % set up the initial array

nspike = size(ISILAT);
min_isi = 2;
if(nspike(2) > 0 && (do_spike == 1))
    spikesumm = spike_sum(DFILE, vsmo, ISILAT, ts1, CONTROL(sf).Rmp, CONTROL(sf).filename, CONTROL(sf).recbeg, CONTROL(sf).spike.deadwin, min_isi, 5, 10, plot_flag);
    % summary spike information:
    if(~isempty(spikesumm))
        n = size(spikesumm.vmax);
        spike.timestamp = datestr(now);
        spike.edit = 0; % not editable...
        spike.n = n(2);
        spike.measthresh = mean(spikesumm.measthresh); % fraction of height from rest to spike peak where measurement was made
        athr = find(~isnan(spikesumm.vthr));
        if(isempty(athr))
            spike.vthresh=NaN;
            spike.vthrmin=NaN;
        else
            spike.vthresh = mean(spikesumm.vthr(athr)); % true threshold (2nd derivative method)
            spike.vthrmin = min(spikesumm.vthr(athr)); % min value (usually the first spike...)
        end;
        spike.deadwin = mean(spikesumm.deadwin);
        spike.minisi = mean(spikesumm.minisi);
        spike.max=mean(spikesumm.vmax); spike.max_sd = std(spikesumm.vmax);
        spike.apht = mean(spikesumm.apht); spike.apht_sd = std(spikesumm.apht);
        spike.half_ht = mean(spikesumm.half_ht); spike.half_ht_sd = std(spikesumm.half_ht);
        spike.min=mean(spikesumm.vmin); spike.min_sd = std(spikesumm.vmin);
        spike.dvmax = mean(spikesumm.dvmax); spike.dvmax_sd = std(spikesumm.dvmax);
        spike.dvmin = mean(spikesumm.dvmin); spike.dvmin_sd = std(spikesumm.dvmin);
        dvvmax = spikesumm.dvvmax(find(~isnan(spikesumm.dvvmax)));
        spike.dvvmax = mean(dvvmax); spike.dvvmax_sd = std(dvvmax);
        dvvmin = spikesumm.dvvmin(find(~isnan(spikesumm.dvvmin)));
        spike.dvvmin = mean(dvvmin); spike.dvvmin_sd = std(dvvmin);
        hwa = spikesumm.hwa(find(~isnan(spikesumm.hwa))); hwb = spikesumm.hwb(find(~isnan(spikesumm.hwb)));
        spike.hwa = mean(hwa); spike.hwa_sd = std(hwa); % eliminate the ones where coulnd't get a measurement
        spike.hwb = mean(hwb); spike.hwb_sd = std(hwb);
        spike.ahpdepth =  mean(spikesumm.ahpdepth); spike.ahpdepth_sd = std(spikesumm.ahpdepth);
        CONTROL(sf).spike=spike; % make it a substructure
        % spike train stuff...
        if(plot_flag >= 0)
            hfig  = newfigure('CC_IV_spiketrain');
            set(hfig, 'units', 'pixels');
            sp = get(hfig, 'Position'); % we need to resize the graph
            set(hfig, 'Position', [sp(1), sp(2), 500 350]);
            hp = multiplot(hfig, 2, 3, 1, [], [8 8], 1);
            set(hfig, 'Renderer', 'Painters'); % vector, best for illustrator transport
            for i = 1:length(hp)
                set(hp(i), 'ytickmode', 'auto');
                set(hp(i), 'yticklabelmode', 'auto');
                set(hp(i), 'xtickmode', 'auto');
                set(hp(i), 'xticklabelmode', 'auto');
            end;

            axes(hp(1));
            plot(spikesumm.lat, spikesumm.dvmax, 'gx');
            xlabel('T (ms)');
            ylabel('dvmax (mV/ms) - rise');
            meandvmax = mean(spikesumm.dvmax(find(spikesumm.lat > 100)));

            axes(hp(3));
            plot(spikesumm.lat, spikesumm.dvmin, 'rx');
            xlabel('T (ms)');
            ylabel('dvmin (mV/ms) - fall');
            meandvmin = mean(spikesumm.dvmin(find(spikesumm.lat > 100)));

            axes(hp(2));
            plot(spikesumm.lat(find(~isnan(spikesumm.hwa))), hwa, 'go');
            xlabel('T (ms)');
            ylabel('hwa (ms) - rising half-width');
            meanhwa = mean(spikesumm.hwa(find(~isnan(spikesumm.hwa(find(spikesumm.lat > 100))))));

            axes(hp(4));
            plot(spikesumm.lat(find(~isnan(spikesumm.hwb))), hwb, 'ro');
            xlabel('T (ms)');
            ylabel('hwb (ms) - falling half-width');
            meanhwb = mean(spikesumm.hwb(find(~isnan(spikesumm.hwb(find(spikesumm.lat > 100))))));

            axes(hp(5));
            axis('off');

            axes(hp(6));
            r_uncomp = CONTROL(sf).access;
            axis([0,1,0,1])
            axis('off')
            text(0,0.80,sprintf('%-12s R[%d:%d] %-8s VCh: %d ICh: %d', ...
                DFILE.filename, DFILE.frec, DFILE.lrec, CONTROL(sf).protocol, vch, ich), 'Fontsize', 9);
            text(0,0.65,sprintf('dVMax: %8.1f dvMin: %8.1f (mV/msec)', meandvmax, meandvmin), 'FontSize', 9);
            text(0,0.50,sprintf('HWA: %8.3f HWB: %8.3f (msec)', meanhwa, meanhwb), 'FontSize', 9);
            text(0,0.35,sprintf('Ihold:%6.2f %s    RMP: %6.2f %s, Rin: %8.3f M\\Omega', ...
                CONTROL(sf).iHold,CONTROL(sf).I_Unit, CONTROL(sf).Rmp, CONTROL(sf).V_Unit, CONTROL(sf).Rin), 'FontSize', 9);
            text(0.75,0.035,sprintf('%s',date),'FontSize',9);
        end;
    end;
end;
CONTROL(sf).spike=spike;



% for measurements from currents and voltages, use data only when spikes DO NOT appear
hypflag = 0;
switch protocol
    case all_hyps
        mean_volt_nspk = mean_volt;
        mean_cur_nspk = mean_cur;
        final_volt_nspk = final_volt;
        peak_volt_nspk = peak_volt;
        peak_tlist_nspk = peak_tlist;
        hypflag = 1;
    otherwise
        nfs=find(isnan(first_spike)); % get the traces with NO spikes.
        mean_volt_nspk = mean_volt(nfs);
        mean_cur_nspk = mean_cur(nfs);
        final_volt_nspk = final_volt(nfs);
        peak_volt_nspk = peak_volt(nfs); % always do "peak" for traces only where NO spikes occur
        peak_tlist_nspk = peak_tlist(nfs);
end

ivresult.timestamp = datestr(now); % place a time stamp on it.
ivresult.edit = 0; % cannot edit - is data calculation result
ivresult.iss = mean_cur; % is mean current
ivresult.vss = mean_volt; % is mean voltage
ivresult.iss_n = mean_cur_nspk; % is mean current excluding spikes
ivresult.vss_n =mean_volt_nspk; % is mean voltage excliding traces with spikes
ivresult.vpk_n = peak_volt_nspk; % is peak voltage excluding traces with spikes (plot against iss_n)
iv.mean=[mean_cur'; mean_volt'];
iv.mean_nspk=[ mean_cur_nspk'; mean_volt_nspk'];
iv.peak_nspk=[ mean_cur_nspk'; peak_volt_nspk'];
iv.mv_fsl = mvoltfsl;
iv.fsl = fsl;
iv.mv_fisi = mvoltfisi;
iv.fisi = fisi;
iv.nr_spk = nr_spikes;

QueMessage('CCIV - Analyzing IV curve (polynomial)');
% make sure cur/volt are "organized..."
[mc,index] = sort(mean_cur);
mv=mean_volt(index);
[mcn, index] = sort(mean_cur_nspk);
pvn=peak_volt_nspk(index);

n_mcn = (mcn-mean(mcn))./std(mcn); % normalize before the fit, then correct later.
pfit = polyfit(n_mcn, pvn, 5); % fit peak voltage with a polynomial line
peak_fit = polyval(pfit, n_mcn);
Rin = diff(peak_fit)./diff(mcn); % use max slope as the input resistance
[C_Rin, Rin_pts] = max(Rin);
if((Rin_pts >= 2) && (Rin_pts < size(Rin, 1)-1))
    CONTROL(sf).Rin = mean(Rin(Rin_pts-1:Rin_pts+1))*1000; % convert from Gohm to Mohm
else
    CONTROL(sf).Rin = Rin(Rin_pts)*1000;
end
%
% do the same for the mean data - used to compute the rectification ratio
%
n_mc = (mc - mean(mc))./std(mc); % normalize prior to doing fit
mfit = polyfit(n_mc, mv, 5);
mean_fit = polyval(mfit, n_mc);
% now extract the Zero Current Potential from this data...
mcd=find(diff(mc)==0);
mc(mcd)=mc(mcd)-0.0001;

CONTROL(sf).Vm0 = polyval(mfit, 0);
%
% now calculate first derivatives from polynomial at 2 symmetrical points
vck1 = 100; % in pA, around rest...
ivslope1 = 5*mfit(1)*vck1^4 + 4*mfit(2)*vck1^3+3*mfit(3)*vck1^2+2*mfit(4)*vck1 + mfit(5);
vck2=-vck1;
ivslope2 = 5*mfit(1)*vck2^4 + 4*mfit(2)*vck2^3+3*mfit(3)*vck2^2+2*mfit(4)*vck2 + mfit(5);

% calculate the mean input resistance in the region from the RMP to spike threshold
ivdepol=NaN;
if(~isnan(CONTROL(sf).spike.vthrmin))
    dep=find(mv>=(CONTROL(sf).Vm0 - 1) & mv < CONTROL(sf).spike.vthrmin); % note the "wiggle" room
    mv_dep=mv(dep);
    mc_dep=mc(dep);
    if(length(mv_dep) > 2) % there must be some room to test this...
        [ivdepol b r p] = linreg(mc_dep, mv_dep); % average is slope if it can be done.
    end;
end;

% store the results - representing slope input resistances at the appropriate levels from rest
ivratio = ivslope2/ivslope1;
iv.timestamp = datestr(now); % analysis time
iv.edit = 0; % not editable - is data calc. result
iv.ivpos = vck1;
iv.ivratio = ivratio;
iv.ivplus = ivslope1;
iv.ivminus = ivslope2;
iv.ivdepol = ivdepol; % slope in region from RMP to just spike threshold
%
% do a curve fit to the more complex equation of a Boltzman with linear term added (and x/y swapped)
% to see if its any different. This is really only valid for type II cells.
% Fit is done in segmented fashion to minimize interactions between leak and v-dep conductance terms.
%
Vr = 0;
T=32.0;	% temp, deg C for typical slice experiment
params = [Vr T];
alpha=(10^-3)*9.648*10^4/(8.3143*(273.16+params(2)));	%1/25mV  =F/RT
maxiter = 150;
if(isempty(CONTROL(sf).Rin))
    % if variables must be empty - Create them, but make them non-existent
    CONTROL(sf).Rin=NaN;
    CONTROL(sf).Rmp=NaN;
end;
do_curr_boltz = 0; % temporarily remove this section.
if(do_curr_boltz)
    QueMessage('CCIV - Fitting IV with Boltzman');
    [imax, jmax] = max(mean_cur);
    gmax = imax/mean_volt(jmax);
    %            gmax      z     Vhalf   gleak                  vleak
    a0 = 0; a1 = gmax; k1 = 4; v1 = -84; x1  = -30;
    [ssiv_par, ssiv_chisq, ssiv_niter, ssiv_fit] = mrqfit('curr_boltzmann', [a0 a1 k1 x1 v1], mean_volt, mean_cur, [], ...
        [0 1 1 1 0], [], [], 50, []);

    iv.gmax = -ssiv_par(2); iv.z = ssiv_par(3); iv.v05 = ssiv_par(5); iv.Gleak = 0; iv.Vleak = 0; iv.fiterr = ssiv_chisq;
    iv.T=T;
    iv.Vr=Vr;
    Vmeas = 5;
    Gplus = 0; b_Gleak = 0;
    %zeRT=alpha*b_z;
    %E = exp(-(mean_volt-b_v05)*zeRT);
    %D1 = (-b_gmax./(1+E)).*(1-(((mean_volt-Vr).*E.*zeRT)./(1+E)))+b_Gleak;
    %sort data before interpolation...
    %[v_sort, index] = sort(mean_volt); % sort by voltage so that ascending order is maintained
    %d1_sort=D1(index);
    %Vm = CONTROL(sf).Rmp+Vmeas;
    %Ep = exp(-(Vm-b_v05)*zeRT);
    %Gplus = (-b_gmax/(1+Ep))*(1-(((Vm-b_v05)*Ep*zeRT)/(1+Ep)))+b_Gleak;
    %Gv = interp1(v_sort, d1_sort, Vm, 'spline');
    %disp(sprintf('Gplus: %8.3f  Gv: %8.3f', Gplus, Gv));
    %Gv
    %b_Gleak
    % save all the parameters here too
    iv.Vmeas = Vmeas;
    iv.Gv = Gplus; % save the conductance measured at that relative voltage
    iv.cratio = Gplus/b_Gleak; % compute the ratio ...
else
    iv.gmax = 0; iv.z = 0; iv.v05 = 0; iv.Gleak = 0; iv.Vleak = 0; iv.fiterr = 0;
    iv.T=0;
    iv.Vr=0;
    iv.Vmeas = 0;
    iv.Gv = 0; % save the conductance measured at that relative voltage
    iv.cratio = 1; % compute the ratio ...
end

% measure the size of the h current - difference between peak and mean voltage (i.e., size
% of sag) at a fixed voltage for the peak current (around -90 mV; use range of +/- 5 mV; i.e., -85 to -95 mV).
%
VSag = -100; % mV
VBand = 10; % mV
vplu = VSag+VBand;
vmin = VSag-VBand;
sag_recs = find(peak_volt >= vmin & peak_volt <= vplu);

if(isempty(sag_recs))
    sag_recs = 1; % pick the first one in the list and try it...
end;
    b = mean(mean_volt(sag_recs));
    a = mean(peak_volt(sag_recs));
    VmSag = b-a;
    VmSagIndex = b/a;
iv.VSag = VSag;
iv.VBand = VBand;
iv.VmSag = VmSag;
iv.VmSagIndex = VmSagIndex; % later we will fit an exponential to that trace

% do alternative measure of the h current - use the mean iv in last 2 msec of current step
% and compute 3 point regressions over appropriate windows
% the windows are : 0 current (from holding) to 2 steps below
% and max step to 2 steps above.
%
stepi = abs(mean(diff(mean_cur))); % average current step size
half_step = stepi/2;
iH = CONTROL(sf).iHold; % get the holding current (again...)
near_rest = find((mean_cur-iH) <= half_step & (mean_cur-iH) >= -(2*stepi+half_step));
mini = min(mean_cur); % find largest step in negative direction (regardless of holding current...)
farfrom_rest = find(mean_cur >= (mini-half_step) & mean_cur <= (mini + (2*stepi+half_step))); % other window

[Rifit, Riint, Ri_r, Ri_p] = linreg(mean_cur(near_rest), mean_volt(near_rest));
[Rpfit, Rpint, Rp_r, Rp_p] = linreg(mean_cur(farfrom_rest), mean_volt(farfrom_rest));
Ri=Rifit; % slope only...
Rp=Rpfit;
rest_i = mean_cur(near_rest);
rest_v = Rifit*rest_i+Riint;
hyp_i = mean_cur(farfrom_rest);
hyp_v = Rpfit*hyp_i+Rpint;

Rhyp = (Rp*Ri)/(Ri-Rp); % compute the parallel equivalent resistance
% disp(sprintf('Ri: %7.4f  Rhyp: %7.4f', Ri, Rp));
iv.Rrest=1000*Ri; % express these in Mohm
iv.Rmaxi=1000*Rp;
iv.Ghyp=1/Rhyp; % express in nS....

% do a time constant fit on the trace at -175 to -225 pA injection...
QueMessage('CCIV - Fitting membrane TC');
iH = CONTROL(sf).iHold;
mincur = min(min(mean_cur-iH)); % get the min current on the traces
minwin_ctr = 0.30*mincur; % take 20% of that value for the selected traces
band = abs(mincur*0.175);
% disp(sprintf('mincur: %f   minwin: %f  - %f', mincur, minwin_ctr-band, minwin_ctr+band))
tfit = find((mean_cur-iH) >= (minwin_ctr - band) & (mean_cur-iH) <= (minwin_ctr + band));
I_avg = mean(mean_cur(tfit)-iH);
g1=zeros(records, 1);
g2=g1;
this_rmp = g2;
for rec=1:records
    g1(rec) = floor(tdel/RATES(rec))-1;
    if(~isempty(deadwin))
        g2(rec) = floor((tdel+deadwin)/RATES(rec));
    else
        g2(rec) = g1(rec);
        deadwin = 0;
    end
end
time_fit=[]; volt_fit=[];
iv.amp0=NaN; iv.amp1=NaN;  iv.tau1=NaN; iv.amp2=NaN; iv.tau2=NaN; iv.tauerr=NaN; % initialize
if(~isempty(tfit) && any(strmatch(protocol, [all_ivs all_hyps {'drugtestiv'}], 'exact'))) % have data to use
    QueMessage(sprintf('CC_IV: Computing taus on %d records', length(tfit)));
    for i=1:length(tfit)
        rec = tfit(i);
        this_rmp(rec) = mean(VOLTAGE(rec,1:trmp(rec)));
        t_list(rec,:) = [g2(rec):peak_tfit(rec)];
        time_fit = time(rec, t_list(rec, :))-tdel-deadwin;
        % do a single exp fit first
        [a0 a1 tau] = expfit(time_fit, vsmo(rec, t_list(rec,:)));
        if(a1 < 0 || tau < 0)
            a1 = 5; % fake it...
            tau = 5;
        end
        [fpar, chisq, niter, volt_fit1(i,:)] = mrqfit('exponential', [a0 a1 tau], time_fit, vsmo(rec, t_list(rec,:)), [], ...
            [], [], [], 50, []);
        if(isempty(volt_fit1))
            volt_fit1(i,:) = a0 + a1*exp(-time_fit/tau);
        end
        fp1(i,:)=fpar;
        cerr1(i)=chisq;
        amp0 = fpar(1); amp1=fpar(2); tau1 = fpar(3);
        [fpar2, chisq2, niter2, volt_fit2(i,:)] = mrqfit('exponential', [amp0 amp1/2 tau1 amp1/2 tau1/8], time_fit, vsmo(rec, t_list(rec,:)), [], ...
            [], [], [], 50, []);
        amp2 = fpar2(1); amp21=fpar2(2); tau21 = fpar2(3); amp22 = fpar2(4); tau22 = fpar2(5);
        time_fit = time_fit+tdel+deadwin;
        if(tau21 > 50 || tau22 > 50)
            chisq2=5*chisq; % make it a worse fit.
        end
        fp2(i,:)=fpar2;
        cerr2(i)=chisq2;
    end % of for loop over records to be fit

    if(mean(cerr2./cerr1) < 0.8 && amp22 > 0)  % require better fit and positive amplitude
        volt_fit= volt_fit2;
        iv.amp0 = mean(fp2(:,1));
        iv.amp1 = mean(fp2(:,2));
        iv.tau1 = mean(fp2(:,3));
        iv.amp2 = mean(fp2(:,4));
        iv.tau2 = mean(fp2(:,5));
        iv.tauerr = mean(cerr2);
        iv.iavg=I_avg;
    else
        volt_fit= volt_fit1;
        iv.amp0 = mean(fp1(:,1));
        iv.amp1 = mean(fp1(:,2));
        iv.tau1 = mean(fp1(:,3));
        iv.tauerr = mean(cerr1);
        iv.iavg = I_avg;
        volt_fit = volt_fit1;
    end
end % end of outer iff..

ih_tau = zeros(length(sag_recs), 1); % all we care about is the tau
t_sag = zeros(length(sag_recs), 2);
sag_volt_fit = cell(length(sag_recs), 1);
sag_time_fit = cell(length(sag_recs), 1);
for i = 1:length(sag_recs)
    rec = sag_recs(i);
    t_sag(i,:) = [peak_v_index(rec),mean_t2(rec)];
    sag_time_fit{i} = time(rec, (t_sag(i,1):t_sag(i,2)));
 % find the min peak voltage during the hyperpolarizing step
    a0 = VOLTAGE(rec, t_sag(i,1));
    a1 = -5;
    tau = 15; % midpoint tau
    % fit the sag of the iH current now:
    [fpar, chisq, niter, sag_volt_fit{i}] = mrqfit('exponential', [a0 a1 tau], sag_time_fit{i}, vsmo(rec, (t_sag(i,1):t_sag(i,2))), [], ...
    [], [], [], 50, []);
    fprintf(1, 'Ih TauDecay: %7.1f\n', fpar(3));
    ih_tau(i) = fpar(3);
end;
iv.VmSagTau = mean(ih_tau);
fprintf(1, 'b/a = %7.3f \n', iv.VmSagIndex);


% do hyp boltzman fit.
if(hypflag)
    QueMessage('CCIV - Fitting FSL with Boltzmann');
    [mvbfsl, ix] = sort(mvoltfsl);
    mvbx=mvbfsl;
    bfsl = fsl(ix); % sort.
    bfslx=bfsl;
    [minfsl, jmin] = min(bfsl); % find max and min fsls.
    [maxfsl, jmax] = max(bfsl);
    jlong = find(bfsl > 2*minfsl); % set up for restriction - jlong indexes all points > 2 x smallest fsl.
    jl2 = jlong(find(jlong > jmin)); % jl2 is all points for v > min fsl, where fsl > 2 x smallest fsl. EXCLUDE these
    jshort = find(bfsl < 0.5*maxfsl); % set up for restriction - jshort indexes all points < 0.5 x largest fsl.
    js2 = jshort(find(jshort < jmax)); % jl2 is all points for v < max fsl, where fsl < 0.5 x largest fsl. EXCLUDE these
    bfsl(jl2) = NaN;
    mvbfsl(jl2) = NaN;
    bfsl(js2) = NaN;
    mvbfsl(js2) = NaN;
    vfsl = intersect(find(mvbx>-80), find(mvbx < -48));
    bfsl = bfslx(vfsl);
    %    mvbfsl
    %    vfsl
    mvbfsl = mvbx(vfsl);
    %   mvbfsl
    %bfsl = bfsl(~isnan(bfsl));
    %mvbfsl = mvbfsl(~isnan(mvbfsl));
    if(length(bfsl) > 3)
        maxiter = 500;
        [fb_a, chisqb_a, niterb_a, fit_ba, errb_a, depb_a] = mrqfit('boltzmann', [0 max(bfsl) -2 -60], mvbfsl, bfsl, [], [1 1 1 1], [], [], maxiter);
        fit_bv = [min(mvbfsl):max(mvbfsl)];
        fit_ba = fb_a(1)+(fb_a(2)-fb_a(1))./(1 + exp(-(- fb_a(4) + fit_bv).* fb_a(3)));

        % store results
        iv.hypb_min = fb_a(1);
        iv.hypb_max = fb_a(2);
        iv.hypb_Vh = fb_a(4);
        iv.hypb_k = fb_a(3);
        iv.hypb_vfit = fit_bv;
        iv.hypb_lfit = fit_ba;
    else
        fit_ba = zeros(1:length(mvbfsl));
        fit_bv = [];
        iv.hypb_min = NaN;
        iv.hypb_max = NaN;
        iv.hypb_Vh = NaN;
        iv.hypb_k = NaN;
        iv.hypb_vfit = NaN;
        iv.hypb_lfit = NaN;
    end;
end;

CONTROL(sf).iv = iv; % make it a substructure - that's how we store results

QueMessage('CCIV - IV analysis complete');

% plot if figure is set
if plot_flag>=0
    fsize = 7;
    msize = 3;
    msg_pos = [0.37 0.00 0.15 0.07];

    hccfig = findobj('Tag', 'CC_IV'); % check for pre-existing window
    if(isempty(hccfig)) % if none, make one
        hccfig = figure('Tag', 'CC_IV', ...
            'Units', 'normalized', ...
            'MenuBar', 'none', ...
            'WindowButtonMotionFcn', 'datac(''mouse_motion'', gcbf);', ...
            'WindowButtonDownFcn', 'datac(''mouse_down'', gcbf);', ...
            'WindowButtonUpFcn', 'datac(''mouse_up'', gcbf);', ...
            'NumberTitle', 'off');
        datac('addwindow', 'CC_IV');
    end
    figure(hccfig); % otherwise, select it
    clf; % always clear the window...
    %Command Menu
    uimenu('Label', 'Close &Window', 'Position', 1, 'callback', 'close(findobj(''Tag'', ''CC_IV''));' );
    uimenu('Label', '&Print', 'Callback', 'print;');
    uimenu('Label', 'Print&Setup', 'Callback', 'printdlg;');
    set(hccfig, 'Name', sprintf('Current Clamp Analysis - File: %s', CONTROL(sf).filename))


    % First plot:
    %plot voltage
    subplot('Position',[0.07,0.52,0.42,0.450]);
    plot(time', VOLTAGE', '-k', ... % data in black, and thin
        time_fit, volt_fit, '-r', ...  % overlay curfit in red
        peak_tlist_nspk, peak_volt_nspk, 'rx'); % peak picked points for analysis with x's
    hold on;
    for i = 1:length(sag_time_fit)
        plot(sag_time_fit{i}, sag_volt_fit{i}, 'b-'); % and blue fit to sag
    end;
    set(gca, 'FontSize', fsize);
    ylabel(CONTROL(sf).V_Unit);
    u=get(gca, 'YLim');
    set(gca,'XTick',[ ]);
    v=get(gca, 'XLim');
    v(2)=tmax;
    set(gca, 'Xlim', v);
    datac_setcrosshair(gca, 'CCIV_info', 'ms', 'mV', msg_pos);

    % title(sprintf('%-12s  R[%i - %i]',DFILE.filename,DFILE.frec, DFILE.lrec));

    %plot currents
    subplot('Position',[0.07,0.4,0.42,0.1]);
    plot(time',CURRENT', '-k');
    set(gca, 'FontSize', fsize);
    v=get(gca, 'XLim');
    v(2)=tmax;
    set(gca, 'Xlim', v);
    ylabel(CONTROL(sf).I_Unit);
    xlabel(CONTROL(sf).T_Unit);
    datac_setcrosshair(gca, 'CCIV_info', 'ms', CONTROL(sf).I_Unit, msg_pos);

    % plot IV
    %        mean_cur_nspk, mean_volt_nspk,'b+', ... % mean voltage - symbol
    %      ssiv_i, ssiv_v, '-g', ... % fit data using boltzman for mean IV
    subplot('Position',[0.07,0.1,0.42,0.25]);
    plot(mean_cur, mean_volt, 'bo', ... % mean voltage - line (or mean_cur_nspk, mean_volt_nspk)
        mean_cur, mean_fit, '-b', ... % fit line using polynomial to the mean IV
        mean_cur_nspk, peak_volt_nspk, 'rx', ... % peak voltage - symbol
        mean_cur_nspk, peak_fit, '-r', ... % peak voltage - line (represents the FIT)
        rest_i, rest_v, '--k',...
        hyp_i, hyp_v, '--k',...
        'MarkerSize', msize);
    set(gca, 'FontSize', fsize);
    xlabel(sprintf('Iinj (%s)',CONTROL(sf).I_Unit));
    ylabel(sprintf('V (%s)',CONTROL(sf).V_Unit));
    h=gca;
    my_legend(h, 'Data-mean', 'bo', 'Fit-Mean', '-b', 'Data-Peak', 'rx', 'Fit-Peak', 'r-');
    set(h, 'FontSize', fsize);
    grid
    % indicate position of Rin measure on IV with an arrow
    text(mean_cur_nspk(Rin_pts), peak_volt_nspk(Rin_pts), '\uparrow', ...
        'HorizontalAlignment', 'Center', ...
        'VerticalAlignment', 'Top');
    datac_setcrosshair(gca, 'CCIV_info', CONTROL(sf).I_Unit, CONTROL(sf).V_Unit, msg_pos);


    % plot spike count in window
    subplot('Position',[0.57,0.72,0.38,0.25]);
    plot(mean_cur,nr_spikes,'-bo', 'MarkerSize', msize);
    grid;
    set(gca, 'FontSize', fsize);
    ylabel(sprintf('Spikes/%4.1f ms', ts2));
    xlabel(CONTROL(sf).I_Unit);
    datac_setcrosshair(gca, 'CCIV_info', CONTROL(sf).I_Unit, 'Sp/Sec', msg_pos);

    % plot first ISI and first spike latency in window
    switch(protocol)
        case {'hyp','hyps', 'ap-hyp', 'ap_hyp'}
            subplot('Position',[0.57,0.4,0.38,0.25]);
            plot(mvoltfsl, fsl, '-rx', mvoltfisi,fisi,'-bo', 'MarkerSize', msize)
            if(hypflag && ~isempty(fit_bv))
                hold on
                plot(fit_bv, fit_ba, '-g');
            end;
            set(gca, 'FontSize', fsize);
            grid;
            ylabel(sprintf('Latency (%s)', CONTROL(sf).T_Unit));
            xlabel(CONTROL(sf).V_Unit)
            legend('1st spike lat', '1st ISI');
            h=gca;
            set(h, 'FontSize', fsize)
            datac_setcrosshair(gca, 'CCIV_info', CONTROL(sf).V_Unit, 'ms', msg_pos);
        otherwise
            subplot('Position',[0.57,0.4,0.38,0.25]);
            plot(mcurfsl, fsl, '-rx', mcurfisi,fisi,'-bo', 'MarkerSize', msize)
            set(gca, 'FontSize', fsize);
            grid;
            ylabel(sprintf('Latency (%s)', CONTROL(sf).T_Unit));
            xlabel(CONTROL(sf).I_Unit);
            if(~isempty(fisi))
                legend({'1st spike lat', '1st ISI'});
            else
                legend('1st spike lat');
            end;
            h=gca;
            set(h, 'FontSize', fsize);
            datac_setcrosshair(gca, 'CCIV_info', CONTROL(sf).I_Unit, 'ms', msg_pos);
    end
    % plot isi vs spike latency
    if(~isempty(ISILAT))
        subplot('Position', [0.57, 0.2, 0.38, 0.15]);
        u = size(ISILAT);
        hold on;
        clist = ['k','r','b', 'g','y', 'm', 'c'];
        clsize = length(clist);
        for i=1:u(2)
            lcolor = sprintf('%c', clist(mod(i, clsize)+1));
            sym = sprintf('-%co', lcolor);
            plot(ISILAT(i).lat(1:length(ISILAT(i).lat)-1), ISILAT(i).isi, sym, 'MarkerSize', msize);
        end

        hold off;
        grid
        set(gca, 'FontSize', fsize);
        ylabel(sprintf('ISI (%s)',CONTROL(sf).T_Unit));
        xlabel(sprintf('Latency (%s)',CONTROL(sf).T_Unit));
        datac_setcrosshair(gca, 'CCIV_info', CONTROL(sf).I_Unit, CONTROL(sf).T_Unit, msg_pos);
    end

    % plot textual information abstracted from analysis...
    subplot('Position',[0.55,0.0,0.45,0.14])
    r_uncomp = CONTROL(sf).access;
    axis([0,1,0,1])
    axis('off')
    text(0,0.8126,sprintf('%-12s R[%d:%d] %-8s VCh: %d ICh: %d', ...
        DFILE.filename, DFILE.frec, DFILE.lrec, CONTROL(sf).protocol, vch, ich), 'Fontsize', 8);
    text(0,0.675,sprintf('Solution:%-12s  gain:%4.1f  LPF:%4.1f kHz', CONTROL(sf).solution, DFILE.igain, DFILE.low_pass(1)), 'FontSize', 7);
    text(0,0.5,sprintf('Ra %5.1f M\\Omega   IVRatio: %6.3f   \\tau1: %6.3f  \\tau2 %6.3f ms', ...
        CONTROL(sf).access, CONTROL(sf).iv.ivratio, CONTROL(sf).iv.tau1, CONTROL(sf).iv.tau2), 'FontSize', 7);
    text(0,0.325,sprintf('Ihold:%6.2f %s    RMP: %6.2f %s, Rin: %8.3f M\\Omega', ...
        CONTROL(sf).iHold,CONTROL(sf).I_Unit, CONTROL(sf).Rmp, CONTROL(sf).V_Unit, CONTROL(sf).Rin), 'FontSize', 7);
    text(0,0.125,sprintf('Mean@ %3.1f-%3.1f ms,   Peak@ %3.1f-%3.1f ms',m1,m2,p1,p2),'FontSize',7);
    text(0.75,0.035,sprintf('%s',date),'FontSize',7);
    orient landscape

    % control printing and closing of window for automatic runs
    % plot_flag = 0 creates plot and leaves it up
    % plot_flag = 1 creates plot and prints to the printer
    % plot_flag = 2 creates plot and prints it and then closes it
    % checking the pflag box in datac main window is plot_flag = 1;
    if (plot_flag > 0)
        print (hccfig, '-dps2');
    end
    if plot_flag == 2
        close
    end
end

